rm(list=ls())
install.packages(c("dplyr", "ggplot2", "ggpubr", "forcats", "broom", "stargazer", "stringr", "scales", "purrr", "tidyverse", "lavaan"))

library(dplyr)      
library(ggplot2)    
library(ggpubr)     
library(forcats)    
library(broom)      
library(stargazer)
library(stringr)    
library(scales)     
library(purrr)      
library(tidyverse)  
library(lavaan)
# Import the CSV file
setwd("C:/Users/bifur/My Drive/Northwestern/My research/Czech Republic/Data")
df <- read_csv("cleaned_df.csv")

# Attention Check Exclusively
df_attf <- df
df<- df %>% filter(attention_check==1)




## Normalize the composite score to a 0 to 1 range and reverse it


df <- df %>%
  rowwise() %>%
  mutate(composite_score = mean(c_across(starts_with("legitimacy_likert")), na.rm = TRUE)) %>%
  ungroup()


df <- df %>%
  mutate(normalized_composite_score = 1 - ((composite_score - min(composite_score, na.rm = TRUE)) /
                                             (max(composite_score, na.rm = TRUE) - min(composite_score, na.rm = TRUE))))



## DESCRIPTIVE STATS ON LEGITIMACY

# Factoring Populism Groups

df$populism_groups <- factor(df$populism_groups, levels = c("Low Populism Score", "Medium Populism Score", "High Populism Score"))



# Grayscale colors for three groups
gray_colors <- c("gray30", "gray60", "gray90")

# Common theme adjustments for all plots
common_theme <- theme(
  plot.title = element_text(size=17, face="italic", hjust=0, margin=margin(t=0, r=0, b=0, l=-10)),  # title with slightly larger size
  axis.text = element_text(size=15),                            # Larger axis text
  axis.title = element_text(size=12),                           # Larger axis title
  legend.text = element_text(size=13),                          # Larger legend text
  legend.title = element_text(size=12),                         # Larger legend title
  plot.margin = unit(c(1, 1, 1, 1), "cm")                       # Increase plot margins
)

# Plot for legitimacy_likert_1
plot_desc_leg1 <- ggplot(df, aes(x=legitimacy_likert_1)) +
  geom_bar(aes(y=(..count..)/sum(..count..), fill=populism_groups)) +
  ylab("Proportion") + xlab("") + labs(fill="") +
  ggtitle("If the Constitutional Court were to issue many decisions\nwith which the majority of the population of the Czech Republic disagrees,\nit would perhaps be better to abolish the Constitutional Court altogether.") +
  scale_fill_manual(values=gray_colors) +
  scale_x_discrete(limits=c("Strongly Disagree", "", "", "", "Strongly Agree")) +
  theme_bw() + common_theme

# Plot for legitimacy_likert_2
plot_desc_leg2 <- ggplot(df, aes(x=legitimacy_likert_2)) +
  geom_bar(aes(y=(..count..)/sum(..count..), fill=populism_groups)) +
  ylab("Proportion") + xlab("") + labs(fill="") +
  ggtitle("The right of the Constitutional Court to decide on some types\nof disputed questions should be limited.") +
  scale_fill_manual(values=gray_colors) +
  scale_x_discrete(limits=c("Strongly Disagree", "", "", "", "Strongly Agree")) +
  theme_bw() + common_theme

# Plot for legitimacy_likert_3
plot_desc_leg3 <- ggplot(df, aes(x=legitimacy_likert_3)) +
  geom_bar(aes(y=(..count..)/sum(..count..), fill=populism_groups)) +
  ylab("Proportion") + xlab("") + labs(fill="") +
  ggtitle("The Constitutional Court can mostly be trusted\nto make decisions that are right for the entire country.") +
  scale_fill_manual(values=gray_colors) +
  scale_x_discrete(limits=c("Strongly Disagree", "", "", "", "Strongly Agree")) +
  theme_bw() + common_theme

# Plot for legitimacy_likert_4
plot_desc_leg4 <- ggplot(df, aes(x=legitimacy_likert_4)) +
  geom_bar(aes(y=(..count..)/sum(..count..), fill=populism_groups)) +
  ylab("Proportion") + xlab("") + labs(fill="") +
  ggtitle("The Constitutional Court meddles too much in politics.") +
  scale_fill_manual(values=gray_colors) +
  scale_x_discrete(limits=c("Strongly Disagree", "", "", "", "Strongly Agree")) +
  theme_bw() + common_theme

# Plot for legitimacy_likert_5
plot_desc_leg5 <- ggplot(df, aes(x=legitimacy_likert_5)) +
  geom_bar(aes(y=(..count..)/sum(..count..), fill=populism_groups)) +
  ylab("Proportion") + xlab("") + labs(fill="") +
  ggtitle("Constitutional Court judges who consistently make decisions\ncontrary to what the majority wishes should be removed from office.") +
  scale_fill_manual(values=gray_colors) +
  scale_x_discrete(limits=c("Strongly Disagree", "", "", "", "Strongly Agree")) +
  theme_bw() + common_theme

# Plot for legitimacy_likert_6
plot_desc_leg6 <- ggplot(df, aes(x=legitimacy_likert_6)) +
  geom_bar(aes(y=(..count..)/sum(..count..), fill=populism_groups)) +
  ylab("Proportion") + xlab("") + labs(fill="") +
  ggtitle("Decisions of the Constitutional Court favor some groups\nmore than others.") +
  scale_fill_manual(values=gray_colors) +
  scale_x_discrete(limits=c("Strongly Disagree", "", "", "", "Strongly Agree")) +
  theme_bw() + common_theme

# Arrange plots in a grid
figure <- ggarrange(plot_desc_leg1, plot_desc_leg2, plot_desc_leg3, plot_desc_leg4,
                    plot_desc_leg5, plot_desc_leg6,
                    ncol = 2, nrow = 3)

# Save the figure with specified size and high resolution
ggsave("populism_legitimacy_measures.png", figure, width = 18, height = 18, units = "in", dpi = 300)



## PARTY LEGITIMACY VISUALIZATION



party_counts <- df %>% 
  count(party_closest) %>% 
  arrange(desc(n))

threshold <- quantile(party_counts$n, probs = 0.1)

df <- df %>%
  mutate(party_closest_to_cats = ifelse(party_closest %in% party_counts$party_closest[party_counts$n <= threshold], 
                                        "Other", 
                                        party_closest))

df <- df %>%
  mutate(party_closest_to_cats = fct_relevel(party_closest_to_cats, "Other"))

df$legitimacy_bin <- cut(df$normalized_composite_score, breaks = seq(0, 1, by = 0.05), include.lowest = TRUE, labels = FALSE)

# Count occurrences for each combination of legitimacy bin and party
df_summary <- df %>%
  group_by(legitimacy_bin, party_closest_to_cats) %>%
  summarise(count = n()) %>%
  ungroup()

# Convert legitimacy_bin to numeric for continuous scale
df_summary$legitimacy_bin <- as.numeric(df_summary$legitimacy_bin)

# Define grayscale color palette
unique_parties <- unique(df$party_closest_to_cats)
num_parties <- length(unique_parties)
gray_colors <- gray.colors(num_parties, start = 0.2, end = 0.8)

# Create the plot using ggplot2
legitimacy_by_party <- ggplot(df_summary, aes(x = legitimacy_bin, y = count, fill = party_closest_to_cats)) +
  geom_bar(stat = "identity", position = "stack") +
  scale_fill_manual(values = gray_colors) +
  scale_x_continuous(breaks = c(1, 6, 11, 16, 20), labels = c("0", "0.25", "0.5", "0.75", "1")) +
  labs(x = "Perceptions of Judicial Legitimacy (Composite)",
       y = "Count",
       fill = "Party Closest To") +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    legend.position = "right"
  ) +
  guides(fill = guide_legend(reverse = TRUE))

ggsave("figure.png", legitimacy_by_party, width = 18, height = 10, units = "in", dpi = 300)



## VISUALIZING PRE-TREATMENT RELATIONSHIPS

# ADD Incumbent Party Variable 

df <- df %>%
  mutate(incumbent_closest = ifelse(party_closest %in% c("SPOLU", "Pirates or Mayors"), 1, 0))

# Define the three models

model_1 <- lm(feeling_cc ~ gender_female + high_school + social_class_normalized, data = df)
model_2 <- lm(feeling_cc ~ gender_female + high_school + social_class_normalized + urban + czech + news_comp + vote_binary, data = df)
model_3 <- lm(feeling_cc ~ gender_female + high_school + social_class_normalized + urban + czech + news_comp + vote_binary+populism_score, data = df)
model_4 <- lm(feeling_cc ~gender_female + high_school + social_class_normalized + urban + czech + news_comp + vote_binary+populism_score+incumbent_closest, data = df)

# Extract coefficients and standard errors from each model
tidy_model_1 <- tidy(model_1, conf.int = TRUE) %>% mutate(model = "Model 1")
tidy_model_2 <- tidy(model_2, conf.int = TRUE) %>% mutate(model = "Model 2")
tidy_model_3 <- tidy(model_3, conf.int = TRUE) %>% mutate(model = "Model 3")
tidy_model_4 <- tidy(model_4, conf.int = TRUE) %>% mutate(model = "Model 4")

# Combine results from all models into one data frame
combined_results <- bind_rows(tidy_model_1, tidy_model_2, tidy_model_3, tidy_model_4)

# Filter out the intercept for plotting
combined_results <- combined_results %>% filter(term != "(Intercept)")

# Recode the term names for better readability
combined_results <- combined_results %>%
  mutate(term = recode(term,
                       "vote_binary" = "Voted (Yes)",
                       "urban" = "Urban",
                       "social_class_normalized" = "Social Class",
                       "populism_score" = "Populist Attitudes",
                       "news_comp" = "General News Exposure",
                       "high_school" = "Education",
                       "gender_female" = "Gender (Female)",
                       "czech" = "Ethnicity (Czech)",
                       "incumbent_closest" = "Support Incumbent"
  ))

# Reverse the levels of 'model'
combined_results$model <- factor(combined_results$model, levels = c("Model 4", "Model 3", "Model 2", "Model 1"))


shape_values <- c("Model 4" = 10, "Model 3" = 15, "Model 2" = 17, "Model 1" = 16)

# Plot with corrected ordering
ggplot(combined_results, aes(x = estimate, y = term, shape = model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray70", size = 0.5) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.3, 
                 position = position_dodge(width = 0.7), color = "gray50") +  
  geom_point(size = 2, position = position_dodge(width = 0.7)) +
  
  labs(x = "Estimated Effect on Feeling Thermometer (Constitutional Court)", y = "", shape = "Model") +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_line(color = "gray90", size = 0.2),
    panel.grid.minor.x = element_line(color = "gray90", size = 0.1),
    axis.title.x = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 12),
    plot.caption = element_text(size = 10),
    plot.background = element_rect(fill = "white", color = NA)
  ) +
  scale_shape_manual(values = shape_values, guide = guide_legend(reverse = TRUE)) +
  scale_x_continuous(limits = c(-0.75, 0.3))



ggsave("descriptive_plot.png", width = 8, height = 6)


## MAIN EFFECTS VISUALIZATION



model <- lm(normalized_composite_score ~ populism + partisanship + irregular + feeling_cc, data = df)
tidy_model <- tidy(model, conf.int = TRUE)

plot_data <- tidy_model %>%
  filter(term %in% c("populism", "partisanship", "irregular"))

# Rename terms for better visualization
plot_data$term <- factor(plot_data$term,
                         levels = c("populism", "partisanship", "irregular"),
                         labels = c("Populism Frame", "Partisanship Frame", "Procedural Regularity Frame"))

ggplot(plot_data, aes(x = estimate, y = term)) +
  geom_point(size = 2, color = "black") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2, color = "gray50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray70", size = 0.5) +
  labs(x = "Estimated Effect on Judicial Legitimacy", y = "") +
  theme_minimal() +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_line(color = "gray90", size = 0.2),  # Lighter and thinner major grid lines
    panel.grid.minor.x = element_line(color = "gray90", size = 0.1),  # Lighter and thinner minor grid lines
    axis.title.y = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)  # White background
  )

ggsave("coefficient_plot.png", width = 8, height = 6)






## VISUALIZING INTERACTION EFFECTS


### Combined Plot

run_interaction_model <- function(data, frame_var, interaction_vars) {
  results <- lapply(interaction_vars, function(interaction_var) {
    # Remove rows with NA values in relevant columns
    complete_data <- data[complete.cases(data[, c("normalized_composite_score", frame_var, interaction_var, "feeling_cc")]), ]
    
    # Construct the formula dynamically
    formula <- as.formula(paste("normalized_composite_score ~", frame_var, "*", interaction_var, "+ feeling_cc"))
    
    # Run the linear model
    model <- lm(formula, data = complete_data)
    
    # Extract the tidy summary with confidence intervals
    tidy(model, conf.int = TRUE) %>%
      filter(str_detect(term, ":")) %>%
      mutate(frame = frame_var, interaction = interaction_var)
  })
  
  bind_rows(results)
}

# Define interaction variables for each frame
interaction_vars_populism <- c("prior_exposure", "gender_female", "social_class", "vote_binary", "newsminuscourt", "populism_score")
interaction_vars_partisanship <- c("prior_exposure", "gender_female", "social_class", "vote_binary", "newsminuscourt", "high_school")
interaction_vars_irregularity <- c("vote_binary", "newsminuscourt", "prior_exposure", "high_school", "gender_female", "social_class")

# Run models and extract results for each frame
results_populism <- run_interaction_model(df, "populism", interaction_vars_populism)
results_partisanship <- run_interaction_model(df, "partisanship", interaction_vars_partisanship)
results_irregularity <- run_interaction_model(df, "irregular", interaction_vars_irregularity)

# Combine results into a single data frame
combined_results <- bind_rows(results_populism, results_partisanship, results_irregularity) %>%
  mutate(frame = factor(frame, 
                        levels = c("populism", "partisanship", "irregular"),
                        labels = c("Populism Frame", "Partisanship Frame", "Procedural Irregularity Frame")),
         interaction = factor(interaction, 
                              levels = c("prior_exposure", "gender_female", "social_class", "vote_binary", "newsminuscourt", "populism_score", "high_school"),
                              labels = c("Prior Exposure to Court News", "Gender (Female)", "Social Class", "Voted (Yes)", "General News Exposure", "Populist Attitudes", "Education")))

# Create the combined interaction effect plot with different shapes for each frame

ggplot(combined_results, aes(x = estimate, y = interaction, shape = frame)) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.1, position = position_dodge(width = 0.5), color = "gray50") +
  geom_point(size = 2, position = position_dodge(width = 0.5), color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray70", size = 0.5) +
  labs(x = "Estimated Interaction Effect on Judicial Legitimacy", y = "", shape = "Frame") +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_line(color = "gray90", size = 0.2),  
    panel.grid.minor.x = element_line(color = "gray90", size = 0.1),  
    axis.title.x = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 12),
    plot.caption = element_text(size = 10),
    plot.background = element_rect(fill = "white", color = NA)  # 
  ) +
  scale_shape_manual(values = c(16, 17, 15)) +
  scale_x_continuous(limits = c(-0.15, 0.15))


ggsave("interactions_combined.png", width = 8, height = 6)


### MAIN EFFECTS TABLE


library(stargazer)

# Assume df is your dataframe
# Run the regressions with normalized_composite_score as the outcome variable

model1a <- lm(normalized_composite_score ~ populism, data = df)
model1b <- lm(normalized_composite_score ~ partisanship, data = df)
model1c <- lm(normalized_composite_score ~ irregular, data = df)

model2a <- lm(normalized_composite_score ~ populism + feeling_cc, data = df)
model2b <- lm(normalized_composite_score ~ partisanship + feeling_cc, data = df)
model2c <- lm(normalized_composite_score ~ irregular + feeling_cc, data = df)






# Create the LaTeX table using stargazer
stargazer(model1a, model1b, model1c, model2a, model2b, model2c, 
          type = "text", #latex 
          title = "Regression Table for Baseline Models",
          dep.var.labels = c("Normalized Composite Score"),
          covariate.labels = c("Populism", "Partisanship", "Irregular", "Feeling CC"),
          omit.stat = c("f", "ser"), # omit F-statistic and Standard Error of the Regression
          no.space = TRUE,
          se = list(summary(model1a)$coefficients[, 2],
                    summary(model1b)$coefficients[, 2],
                    summary(model1c)$coefficients[, 2],
                    summary(model2a)$coefficients[, 2],
                    summary(model2b)$coefficients[, 2],
                    summary(model2c)$coefficients[, 2]))

### INTERACTIONS TABLES


# Create interaction models for Populism
model4a <- lm(normalized_composite_score ~ populism * prior_exposure + feeling_cc, data = df)
model4b <- lm(normalized_composite_score ~ populism * gender_female + feeling_cc, data = df)
model4c <- lm(normalized_composite_score ~ populism * social_class + feeling_cc, data = df)
model4d <- lm(normalized_composite_score ~ populism * vote_binary + feeling_cc, data = df)
model4e <- lm(normalized_composite_score ~ populism * newsminuscourt + feeling_cc, data = df)
model4f <- lm(normalized_composite_score ~ populism * populism_score + feeling_cc, data = df)
model4g <- lm(normalized_composite_score ~ populism * education + feeling_cc, data = df)

# Create interaction models for Partisanship
model5a <- lm(normalized_composite_score ~ partisanship * prior_exposure + feeling_cc, data = df)
model5b <- lm(normalized_composite_score ~ partisanship * gender_female + feeling_cc, data = df)
model5c <- lm(normalized_composite_score ~ partisanship * social_class + feeling_cc, data = df)
model5d <- lm(normalized_composite_score ~ partisanship * vote_binary + feeling_cc, data = df)
model5e <- lm(normalized_composite_score ~ partisanship * newsminuscourt + feeling_cc, data = df)
model5f <- lm(normalized_composite_score ~ partisanship * high_school + feeling_cc, data = df)

# Create interaction models for Procedural Irregularity
model6a <- lm(normalized_composite_score ~ irregular * prior_exposure + feeling_cc, data = df)
model6b <- lm(normalized_composite_score ~ irregular * gender_female + feeling_cc, data = df)
model6c <- lm(normalized_composite_score ~ irregular * social_class + feeling_cc, data = df)
model6d <- lm(normalized_composite_score ~ irregular * vote_binary + feeling_cc, data = df)
model6e <- lm(normalized_composite_score ~ irregular * newsminuscourt + feeling_cc, data = df)
model6f <- lm(normalized_composite_score ~ irregular * high_school + feeling_cc, data = df)


## Populism Table

stargazer(model4a, model4b, model4c, model4d, model4e, model4f, model4g, 
          type = "text", #latex
          title = "Populism Interaction Models",
          dep.var.labels = "Legitimacy Attitudes",
          no.space = TRUE)

## Partisanship Table

stargazer(model5a, model5b, model5c, model5d, model5e, model5f, 
          type = "text", #latex
          title = "Partisanship Interaction Models",
          dep.var.labels = "Legitimacy Attitudes",
          no.space = TRUE)

## Irregularity Table

stargazer(model6a, model6b, model6c, model6d, model6e, model6f, 
          type = "text", #latex
          title = "Procedural Irregularity Interaction Models",
          dep.var.labels = "Legitimacy Attitudes",
          no.space = TRUE)

## Alt. SEM Approach



model <- '
  lativ =~ legitimacy_likert_1 + legitimacy_likert_2 + 
                               legitimacy_likert_3 + legitimacy_likert_4 + 
                               legitimacy_likert_5 + legitimacy_likert_6
  
  lativ ~ populism + partisanship + irregular + feeling_cc
'

fit <- sem(model, data = df)

summary(fit, fit.measures = TRUE)